function [f, sef] = gammafunc_se(T,beta, M, pxest)

theta1 = ones(M,1);
dgammadrho = zeros(M,1);
for m=1:M
    
B = 1;
for t=2:T
temp = beta^(t-1);
temp2 = (beta^(t-1))*( pxest(2,m)^(t-1));
temp3 = (beta^(t-1))*( (pxest(2,m)^(t-2))*(t-1));
B = B+temp;
theta1(m) = theta1(m) + temp2;
dgammadrho(m) = dgammadrho(m) + temp3;
end
end
theta1 = theta1/B;
dgammadrho = dgammadrho/B;

f = theta1;
sef = dgammadrho;

end


